Membrane-Protein Interactions in a Generic Coarse-Grained Model for Lipid Bilayers 
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We study membrane-protein interactions and membrane-mediated protein-protein interactions 
by Monte Carlo simulations of a generic coarse-grained model for lipid bilayers with cylindrical 
hydrophobic inclusions. The strength of the hydrophobic force and the hydrophobic thickness of 
the proteins are systematically varied. The results are compared with analytical predictions of two 
popular analytical theories: The Landau-de Gennes theory and the elastic theory. The elastic theory 
provides an excellent description of the fluctuation spectra of pure membranes and successfully 
reproduces the deformation profiles of membranes around single proteins. However, its prediction 
for the potential of mean force between proteins is not compatible with the simulation data for large 
distances. The simulations show that the lipid-mediated interactions are governed by five competing 
factors: Direct interactions, lipid-induced depletion interactions, lipid bridging, lipid packing, and a 
smooth long-range contribution. The mechanisms leading to "hydrophobic mismatch" interactions 
are critically analyzed. 
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I. INTRODUCTION 

Membrane proteins are integral components of 
biomembranes and account for most of the biological pro- 
cesses that take place in and at membranes Their 
activity often depends on their distribution within the 
membrane 0] • The latter is determined by direct interac- 
tions, but also to a significant extent by indirect interac- 
tions, such as those mediated by the lipid bilayer matrix. 
Lipid-protein interactions are believed to play an impor- 
tant role, e.g., in controlling the aggregation and activ- 
ity of gramicidin channels [|[ and rhodopsin There- 
fore, lipid-mediated interactions between membrane pro- 
teins or more generally membrane inclusions have been 
studied intensely for many decades [EBB]- Natural 
biomembranes are of course complex multicomponent 
systems, and membrane heterogeneities contribute crit- 
ically to the lipid-protein interactions @. However, a 
number of mechanisms have been identified that gener- 
ate lipid-mediated protein interactions already in pure, 
one-component lipid bilayers. 

(i) The mechanism which has been pointed out first 
in the literature is the hydrophobic mismatch in- 
teraction [H, [l(| H3 • It arises in situations where 
the hydrophobic thickness of transmembrane pro- 
teins does not match the equilibrium bilayer thick- 
ness. The proteins then locally compress or ex- 
pand the membrane, and the associated free energy 
penalty depends on their distribution in the mem- 
brane. This may induce protein clustering. The 
effect has been verified experimentally with system- 
atic studies of gramicidin [12| and synthetic model 
peptides [HI, Theoretically, it has been ex- 

plained using different continuum theories for bi- 



layers ha ha p. |ia ha |2fl |2il I22, y , y, |25|, 

HE USUI, Ha Ho, \M Wk- However, it does not 
seem strong enough to fully account for the exper- 
imentally observed clustering of proteins in mem- 
branes [33l |34|. 

(ii) Even in the absence of hydrophobic mismatch, in- 
clusions locally disturb the translational and con- 
formational degrees of freedom of lipids [{J Ef| [H, 
[37l [H, HtJ E| , which leads to local packing in- 
teractions [42|. They have been analyzed theoreti- 
cally by sophisticated mean-field studies of effective 
interactions between fully repulsive inclusions in- 
serted in membranes of fixed thickness [13, [H, Hil . 
These calculations generally predict attractive in- 
teractions at short distances and repulsive interac- 
tions at larger distances. 

(hi) A third class of interactions discussed in the litera- 
ture are fluctuation-induced interactions. Proteins 
locally affect the elastic properties of the mem- 
branes - the bending rigidity and/or the spon- 
taneous curvature - thereby changing the fluc- 
tuation spectrum of the membrane. The corre- 
sponding entropy change depends on the distri- 
bution of the proteins, which leads to Casimir 
forces [H BEE [3. 

(iv) In addition to these general interaction mecha- 
nisms, a number of more specific effects have been 
studied. For example, lipid-mediated interactions 
are observed between proteins that locally induce 
a strong membrane curvature. The sign of these 
interactions is not clear - whereas elastic theories 
predict repulsion, curvature-inducing proteins have 
in fact been found to aggregate in coarse-grained 
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simulations [47|. This indicates again that elastic 
theories alone cannot fully account for the effec- 
tive interactions in such a system. Complex in- 
teraction mechanisms have also been predicted for 
membranes in low-temperature ordered, e.g., tilted 
states [H|. 

In the present article, we study the interplay of differ- 
ent lipid-mcdiated interaction mechanisms between sim- 
ple cylindrical inclusions in one-component lipid bilay- 
ers by computer simulations of a generic coarse-grained 
membrane model. Despite being very simple, our model 
reproduces the main phases and conformationally driven 
phase transitions of pure lipid layers, including struc- 
tures as complex as asymmetric and symmetric rippled 
states [iil . It seems therefore suited to study generic phe- 
nomena that depend on lipid conformations. Specifically, 
we focus on the first two interaction mechanisms listed 
above, the hydrophobic mismatch interaction (i), and the 
packing interactions (ii). The diameters of our inclusions 
are too small to generate measurable Casimir forces, they 
correspond roughly to those of simple /3-helices. The in- 
clusions do not induce curvature, and the bilayer is in the 
fluid state, hence additional interactions (iv) also do not 
contribute. 

In the past decades, a number of computer simula- 
tion studies have focused on protein-lipid interactions, 
using both atomistic and coarse-grained models (see, 
e.g., Refs. [Hoi . l5lj for recent overviews). Simulation 
studies of membrane-mediated protein-protein interac- 
tions are less abundant. Sintes and Baumgaxtner [52| 
have been the first to study lipid-mediated interactions 
in a coarse-grained molecular model. They considered 
purely repulsive cylinders immersed in a bilayer of lipids 
which are head-grafted to opposing flat surfaces. In 
qualitative agreement with the mean-field theories cited 
above [13, 0, [39| , they find an attractive depletion in- 
teraction at close distances followed by a repulsive well. 
Smeijers et al. [53| have studied the aggregation of pro- 
teins with different shapes in fusing vesicles. Very re- 
cently, de Meyer, Venturoli and Smit [5l[ have presented 
a systematic study of lipid-mediated protein interactions 
for varying hydrophobic mismatch in a coarse-grained 
membrane model with soft DPD (dissipative particle) 
interactions. Based on their data, they argue that an 
important factor driving the hydrophobic mismatch in- 
teraction is "hydrophilic shielding", i.e., the influence of 
mismatch on the local arrangement of head groups rela- 
tive to tails. 

The present study is in many resp ect complementary 
to the work of de Meyer et al. [51j. First, our model 
is very different: On the one hand, our coarse-grained 
"lipid" structure is much simpler, on the other hand, we 
use hard core, Lennard-Jones type interactions, similar to 
those used in atomistic or systematically coarse-grained 
models [54[ • This allows us to study the influence of lo- 
cal lipid packing phenomena, which are almost absent 
in systems with soft DPD potentials, and to diagnose 
new factors that might contribute to the hydrophobic 



mismatch interaction, such as local chain ordering. Sec- 
ond, we systematically vary the hydrophobicity of the 
protein and study its influence on the protein-protein in- 
teractions. Third, we relate our simulation results to two 
popular analytical theories of lipid-induced interactions: 
The Landau-de Gennes theory flElfl^. 17|and the elastic 



theory for coupled monolayers (1171217130, EL H- The 
elastic theory turns out to provide an excellent descrip- 
tion of the peristaltic and bending fluctuations of pure 
membranes, and of thickness deformation profiles around 
single proteins. This allows us to analyze the elastic con- 
tribution to the protein-membrane interactions - among 
other things, we identify a mechanism how they may be 
affected by hydrophilic shielding. The simulation results 
for the membrane-mediated protein- protein interactions, 
however, are not compatible with the predictions of the 
elastic theory, especially for large protein distances. In 
this case the simpler Landau-de Gennes theory seems to 
perform better. 

Our paper is structured as follows: After introducing 
the model and the simulation method and briefly rec- 
ollecting the main assumptions and predictions of the 
theories in the next section, we present and discuss our 
simulation results in the section Three. We conclude with 
a brief summary. 



II. MODELS AND METHODS 

A. Simulation Model and Methods 

We employ a simple generic lipid model [55| that 
has been shown to reproduce the characteristic high- 
temperature phase transitions of monolayers [H, H3] and 
bilayers (49[. The lipids are represented by chains of 
seven beads with one head bead of diameter ah followed 
by six tail beads of diameter a t . Beads, that are not di- 
rect neighbors in a chain interact, via a truncated and 
lifted Lennard-Jones potential: 



V LJ (r/a)-V LJ (r c /a) 




if r < r c 
otherwise 
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where a is the mean diameter of the two interacting 
beads, cry = (<7j + <Xj)/2 (i,j = h or t). Head-head and 
head-tail interactions are purely repulsive (r c = a) while 
tail-tail interactions also have an attractive contribution 
(r c = 2a). Within a "lipid" chain beads are connected 
to each other by FENE (Finitely Extensible Nonlinear 
Elastic) springs with the spring potential 
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where tq is the equilibrium distance, Ar max the maximal 
deviation, and cfene the FENE spring constant. In ad- 
dition, the chains are given bending stiffness by means of 
a bond-angle potential 



Vba(0) = £ba(1 - cos(0)). 



(4) 



The aqueous environment of the membrane is modeled 
with "phantom" solvent beads [581 ] . which interact with 
lipids like head beads (a s = ah), but have no interac- 
tions with each other. Much like an implicit solvent, the 
phantom solvent is structureless and does not impart un- 
wanted correlations onto the bilayers, and it is very cheap 
from a computational point of view. At sufficiently low 
temperatures it forces the lipids to spontaneously self- 
assemble into stable bilayers |58j |. 

Specifically our model parameters are [EH, [E3] crh = 
l.lCTt, r = 0.7t7 t , Ar max = 0.2<r t , e FENE = lOOe/of , and 
e ba = 4.7e. The simulations were carried out at con- 
stant pressure P = 2.e/af and temperature fcgT = 1.3e, 
which is well in the fluid phase region of the bilayer. 
(The bilayer undergoes a main transition to a tilted 
gel phase Lpi via a ripple phase Ppi at the tempera- 
ture ksT = 1.2e |4?|.) Comparing the monolayer thick- 
ness, to ~ 3ct(, and the area per lipid, ao ~ 1.36er 2 , 
with the corresponding numbers for real lipid bilayers 
in the fluid L a phase, we find that the values in our 
model roughly reproduce those of DPPC (dipalmitoyl 
phosphatidylcholine) bilayers, if we set a t ~ 6A [59j |. 
By matching the temperatures of the main transition 
(T m = 42°C in DPPC) we can also identify an energy 
scale: e ~ 0.36 • lO" 20 .!. 

The proteins are modeled as cylinders with the radius 
R = 1.5ct(, which have fixed orientation along the bi- 
layer normal (the z-axis). We impose the orientation in 
order to realize the situation considered in the elastic 
theory as closely as possible - proteins cannot respond 
to hydrophobic mismatch by tilting. In a real membrane, 
this would correspond to a situation where the orienta- 
tion of the transmembrane domain of a protein is kept 
fixed by external factors, e.g., geometrical constraints in 
the extramembrane domain. For comparison, selected 
simulations were also conducted for proteins which were 
allowed to tilt (unconstrained orientations). As in other 
model studies |3ll . l5l| , the effect of orientation fluctua- 
tions on the results was found to be rather small (see 
section Three). 

The interaction between proteins and lipid or solvent 
beads has a purely repulsive contribution, which is de- 
scribed by a radially shifted and truncated Lennard- 
Jones potential 



VU 1 



VLj(1) 



if r — <7q < a 
otherwise 



(5) 



where r = y/ x 2 + y 2 denotes the distance of the interact- 
ing partners in the (a;y)-plane, the parameter a is given 
by a = ((T( + <t,)/2 for interactions with beads of type 
i (i = h, t, and s for head, tail, and solvent beads, re- 
spectively), ctq = at, and Vlj has been defined above 



(Eq. ([2])). The direct protein-protein interactions have 
the same potential (0 with a = a t and <j = 2<r t . 

In addition, protein cylinders attract tail beads on a 
"hydrophobic" section of length L. This is described by 
an additional attractive potential that depends on the 
z-distance between the tail bead and the protein center. 
The total potential reads 

Vpt(r, z) = e pt (V rcp (r) + Kttr(r) • Wp(z)) (6) 

with the attractive Lennard- Jones contribution 



Kttr(r) 



- F LJ (2) if r-a < a 
VLj(^) - V LJ (2) if a < r - a < 2a 
otherwise 



(7) 

and a weight function Wp (z) , which is unity on a stretch 
of length L — 2a t and crosses smoothly over to zero over 
a distance cr t at both sides. Specifically, we use 

1 if |;?| < J 

W P (z) = { cos 2 [3/2(|z| -I)} if I < \z\ < 1 + t:/3 (8) 

otherwise 

with / = L/2 — a t . The parameter e p t tunes the strength 
of the lipid-protein interaction, i.e., the hydrophobicity 
of the protein. It was varied between e pt = 1 and e pt = 6. 
The hydrophobicity e p t — 1 was sufficient to trap the cen- 
ter of the protein inside the membrane: The fluctuations 
of its z-position relative to the height of the membrane 
h were of the order ((z protein — h) 2 ) ~ 0.5cr t for all values 
of e p t . The proteins with unconstrained orientation were 
modeled as spherocylinders of length L with essentially 
the same interaction potentials, except that the z-axis is 
replaced by the protein axis, and the variable r by the 
closest distance to the protein. 

The system is studied using Monte Carlo simulations at 
constant pressure and temperature with periodic bound- 
ary conditions in a simulation box of variable size and 
shape: The simulation box is a parallelepiped spanned by 
the vectors (L x , 0, 0), (s yx L x , L y , 0), (s zx L x , s zy L y , L z ), 
and all Li and Sj are allowed to fluctuate. This ensures 
that the membranes have no interfacial tension. The sys- 
tem sizes ranged from 200 to 3200 lipids, and the simula- 
tions were run up to 8 million Monte Carlo steps, where 
one Monte Carlo step corresponds to one Monte Carlo 
move per bead. Moves that alter the simulation box were 
attempted every 50th Monte Carlo step. To generate the 
initial configurations, we set up a perfectly ordered bi- 
layer in the (a;y)-plane with straight chains pointing in 
the z-direction and simulated it until it was equilibrated, 
i.e., all observables of the system fluctuated about the 
equilibrium value and none of the observables shows a 
trend. Typical equilibration times where 1 million Monte 
Carlo steps (4 million Monte Carlo steps in simulations 
where we looked at large-scale height fluctuations) . Due 
to this procedure, the bilayers were oriented in the (xy)- 
plane. 
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FIG. 1: Cross-section snapshot of a model membrane with 
two inclusions of hydrophobic thickness L = 6 and different 
hydrophobicity parameter e pt = 3 (top) and e pt — 6 (bottom). 
Light grey lines show tail bonds, dark circles the heads (not 
to scale), dark cylinders the inclusions. 



deformation profile 



0Ldo(r) = t R 



Ko(R/0' 



(10) 



where R is the radius of the inclusion, £ = yjc/a the 
correlation length, and Kq is the modified Bessel func- 
tion of the second kind. For r£ > 1, Eq. (fTOf can be 
approximated by an exponential decay 



(r) w t R exp 




(11) 



Such exponential laws have often been used to fit data 
from simulations or molecular theories for membranes 
with inclusions [H, H3, IM HI 13] . Eq. ([9]) can also be 
used to deduce an equation for the monolayer thickness 
fluctuation spectrum 



(\^a(q)\ 2 ) = 



k B T 



4(a + cq 2 ) ' 



(12) 



Fig. Q] shows two snapshots of a system containing two 
inclusions whose thickness roughly matches that of the 
bilayer [L = 6a t), with different hydrophobicity param- 
eters e p t- They illustrate the existence of membrane- 
mediated attractive interactions even in the absence of 
hydrophobic mismatch. At low values of e p t, the pro- 
teins touch. At higher values of e p t, they are separated 
by a single lipid layer. 

Before proceeding to a more quantitative discussion of 
the simulation results, we shall now briefly recollect the 
main assumptions and results of the analytical theories 
which we use to analyze our data. 



B. Landau-de Gennes Theory 

One of the oldest theoretical approaches to studying 
lipid-mediated interactions between inclusions is based 
on a Landau-de Gennes expansion of the free energy in 
powers of the li pid area or membrane thickness varia- 
tions fl5l. H5 , [l7l. [60| . In the simplest case, this expansion 
reads [13] 

F, dG = J d 2 r{|(20 LdG ) 2 + ^(2V</> LdG ) 2 } (9) 

with the boundary condition LdG = ifi at the surface 
of the inclusion, where LdG denotes the local deviation 
of the monolayer thickness from its equilibrium value to 
in the unperturbed membrane, and 2(tn + to) is the hy- 
drophobic thickness of the inclusion. The first term in 
Eq. © accounts for the area compressibility fc^ of the 
bilayer Qca = a(2io) 2 ) 7 and the second term penalizes 
spatial thickness variations, i.e., the variable c is taken 
to be positive. Minimizing this free energy for a mem- 
brane containing a single inclusion at r = yields the 



C. Elastic Theory 

Another popular approach is the elastic th eory of cou- 
pled monolayers [H M, US HE H3, SI, M, M, M, H3] • 
Here the membrane is treated as a system of two coupled 
elastic monolayer sheets. The basic structure of the dif- 
ferent theories is very similar - they differ mainly in the 
choice of the boundary conditions and the number of elas- 
tic terms that they include. For example, early theories 
have often disregarded the possibility that the individual 
monolayers may have a spontaneous curvature, whereas 
this is usually accounted for in more recent work. Here 
we shall use a recent version of the elastic theory devel- 
oped by Brannigan and Brown [3lL |32| . which is fairly 
"complete" in the sense that it includes all known elastic 
terms. 

We consider a flat lipid bilayer in the (xy)-p\ane. The 
two constituting monolayers are described by four in- 
dependent fluctuating fields - two accounting for meso- 
scopic bending deformations, and two for the microscopic 
protrusions. We assume that the volume of lipids is lo- 
cally conserved, and that the mesoscopic bilayer height 
and thickness fluctuations and the protrusions basically 
decouple. The latter requirement may seem all-too rigid 
and is actually not imposed in the original model [311 ] , but 
it leads to a considerable simplification of the resulting 
theory and will be justified a posteriori by the fact that 
it describes the fluctuation spectra of our model bilayers 
in an excellent way. 

The quantities of interest for us are the local bilayer 
height h, the local monolayer thickness t, and the local 
field </> cl which denotes just the contribution of bending 
deformations to the thickness, without the microscopic 
protrusions. This field is defined in analogy to the field 
(/> LdG in Eq. (0). With the assumptions mentioned above, 
the spectra of the bilayer height and monolayer thickness 
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fluctuations in Fourier space, (\h(q)\ 2 ) 
given by [3 11 ] 



and (\t(q)\ 2 ), are at infinity, where V 2 = (l/r)d r rd r and 



and 



(l%)| 2 ) 
(\t(q)\ 2 ) 



k c q 4 + 2(k x +^q 2 ) : 
k^T 

k^-Akc^/to + kA/tl 

k B T 
"2(fc A + 7 A? 2 )' 



(13) 
(14) 



where k c and kA are the bending and the compressibility 
modulus of the bilayer, £ is related to the spontaneous 
curvature 64], to is the mean monolayer thickness, and 
the parameters 7a and k\ characterize the protrusions. 
Fitting our simulation data to these theoretical spectra 
allows us to (i) test the validity of the theory and (ii) ex- 
tract the elastic parameters k c , /ca, and £, which we need 
for the subsequent analysis of protein-lipid interactions. 

Within the decoupling approximation, the protrusions 
and the height fluctuations do not contribute to the 
protein- induced membrane deformations. The free en- 
ergy of monolayer thickness deformations can be ex- 
pressed as [HJ 



F, 



el.O 



to 



y(V 2 OI ) 2 + fc G det(c^0 ol ) 



(15) 



where (</> cl + to) is the locally smoothed monolayer thick- 
ness (without the protrusions), and we have introduced 
the spontaneous curvature of the monolayers cq and the 
Gaussian rigidity kc- According to the Gauss-Bonnet 
theorem, the latter only contributes an uninteresting con- 
stant in homogeneous planar sheets and is thus often 
omitted; in the presence of inclusions (holes), however, 
it has to be taken into account [32[ . 

We consider inclusions with radius R that enforce a 
certain membrane thickness 2tn at their surface. To cal- 
culate the deformation profile of the bilayer in the vicinity 
of such an inclusion centered at r = 0, we minimize the 
free energy F e ; i0 with respect to the profile </> cl (r) while 
keeping the membrane deformation at the surface of the 
inclusion fixed, ^surface = j. R rp n j g Y ea( \ s to the Eulcr- 

Lagrange equation 



kA 
k t 2 ' 



to 



+ V 4 



= 



with the boundary conditions 



MR) 
v 2 . 



t R 



kc ,, f , Jr 

R = ^fi ~ 2 ( c + C — 



k c R 



at the surface of the inclusion, and 



<9r<^=l(r)|r 



V^ el (r)| r _ 



(16) 

(17) 
(18) 

(19) 



we have defined t' R = d r <fi al \R. For a single inclusion, 
these equations can be solved analytically, giving [3(3| 



ol (r) = Ai Jo{a+r) + A 2 Y (a+r) 
+ A 3 J {a-r) + A 4 Y (a-r) 



(20) 



with |6E 



a± 



\ 




j l — 1771' (21) 



where Jo{x) and Yq(x) are the zeroth order Bessel func- 
tions of the first and second kind, and the coefficients Ai 
are determined by the boundary conditions. 

The elastic model presented so far only uses material 
properties of free, bulk membranes. Inclusions may lo- 
cally alter the lipid properties, e.g., the lipid volume, the 
lipid ordering etc., which may in turn affect the elastic 
properties (e.g., the equilibrium thickness, the sponta- 
neous curvature, the bending rigidity, the compression 
modulus) of the membrane. In Ref. [32j . Brannigan and 
Brown have demonstrated for the case of varying lipid 
volume that such effects can be incorporated into the 
theory in a relatively straightforward way [32j. 

Here we will discuss this from a more general perspec- 
tive: Consider some scalar quantity q(r) that is distorted 
from its bulk value go by the inclusion and locally alters 
the membrane properties. By symmetry, it will introduce 
two new terms in Eq. (|15[) . 



F, 



1.0 



F n 



with 



= ; ! i\\ —< 

qo 



9o 



(22) 



Here Sq/qo denotes the relative deviation of q, and terms 
that do not depend on (f> cl or that are of higher than 
quadratic order in the deviations <p cl and Sq/qo have been 
disregarded. In general, the additional free energy contri- 
bution F q will change the Euler-Lagrange equations, and 
the solution ([2H)l is no longer valid close to the inclusion. 
The situation however simplifies considerably if we make 
the reasonable assumption that the local distortion 5q(r) 
decays to zero on a length scale which is much smaller 
than the characteristic length scales of the elastic profile. 
For an inclusion centered at r = 0, we can then replace 
cl (r) by tn + t' R (r — R) in Eq. (f22|) . and the free energy 
F q turns into a surface term, 



Fq = t R K X 



2ndr r 



Sq(r) 



+ t' t 



R 90 

, 30 dr2^M!l(X 1 r(r - R) + K 2 ) 

r qo 



(23) 



which only changes the boundary condition (fT8|) : The 
local distortion 5q(r) renormalizes the spontaneous cur- 
vature term in Eq. (TT8"]) according to 



1 



CO = c 



2k r R 



dr Mil (#i r (r _ R ) + K 2 ). (24) 
9o 
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Inserting that into Eq. (|22|) with ([15]) and exploiting the 
Euler Lagrange equation (HHJ), we find that the free en- 
ergy of the deformation is given by 

F ehq = 7rk c R(t R Vl<f> el \ R - 2t' R (ca - ( t R /t )) + const., 



100 



where the constant does not depend on the deformation 
profile. 

This result can readily be generalized to situations 
where several scalar quantities q a (r) affect the membrane 
simultaneously. Within our linear approximation, each of 
them will contribute a separate surface term F qa of the 
form (|23|) and the effect on the renormalized curvature 
term (|24[) will be additive. For given renormalized cur- 
vature Co, Eq. (|25|) still remains valid. 

Our findings agree qualitatively with those of Branni- 
gan and Brown in Ref. [32| , who also conclude that lipid 
volume deviations v(r)/vo at the surface of the inclusion 
effectively renormalize the spontaneous curvature term. 
The actual expression for c~o given in that paper is differ- 
ent from ours. The discrepancy is due to an error in the 
original analysis. 



D. Other Theories 
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FIG. 2: Fourier spectra of height (closed circles) and thickness 
fluctuations (open squares). The dashed line shows a fit of the 
data to the elastic theory |3l[ (Eqs. (fl3|l and CI}). The inset 
shows the thickness data alone with a fit to the Landau-de 
Gennes theory (Eq. <JT2j) ) . 



The elastic theory sketched in the previous section is 
an effective interface theory, i.e., the main degrees of 
freedom are the positions of the fluctuating interfaces 
between the monolayers and the surrounding solvent. 
A number of authors have put forward elastic theories 
which include the local tilt of chains as supplementary 
degrees of freedom in the spirit of the Landau-de Gennes 
theory for smectic liquid crystals [12, 51, [H, [63] • These 
theories introduce new elastic parameters which describe, 
e.g., the splay, twist, and bend modes of the tilt order pa- 
rameter, and their coupling to the interfacial degrees of 
freedom. Due to the difficulty of accessing the additional 
parameters, they are not included in the present analysis. 

Another entirely different type of approach is pursued 
by the molecular theories [H [H, [33, H [H, ||] . They 
account for the chain character of lipids explicitly and 
calculate the packing interactions between proteins and 
membranes by different sophisticated mean field meth- 
ods. Since the results depend on the chain model, which 
is different from ours, they can only be compared to our 
simulation data at a qualitative level (see next section). 

III. SIMULATION RESULTS 

We turn to discuss the simulation results. First, we 
consider the properties of pure bilayers, with no inclu- 
sions. The results confirm the validity of the elastic the- 
ory for our system, and provide values for the elastic 
parameters that can be used subsequently. In the second 
step, we investigate the deformation profiles of a mem- 
brane in the vicinity of single inclusions. Last, we study 



the effective inclusion-inclusion interactions for inclusions 
with different hydrophobic lengths and hydrophobicity 
parameters. 



A. Elastic Properties of the Pure Membrane 

One of the most powerful approaches to looking at 
elastic membra ne prop erties is the analysis of fluctua- 
tion spectra [1^, [7(J O, [zl] • To determine the spectra for 
our membranes, we have basically followed the procedure 
outlined in Ref. [l3] . We have carried out simulations of a 
bilayer containing 3200 lipids (and 24615 solvent beads). 
The system was divided in N x x N y bins in the (xy)-plane 
with N x — N y — 20. In each bin, the z-value of the mean 
head position in the z-direction was determined for both 
monolayers separately. The average and the difference of 
the two values give the bilayer height spectrum h(x,y) 
and monolayer thickness spectrum t(x,y), respectively. 
The two spectra were then Fourier transformed accord- 
ing to 

f** = w£Y,K x >V>~ H{ * mW+9 ' V) > (26) 

x y x,y 

and the values of \h q \ 2 and \t q \ 2 were collected in g-bins 
of size 0.1. The binning was made necessary by the fact 
that the dimensions of the simulation box and hence the 
q- values fluctuate. In each bin, the averages (\h q \ 2 ) and 
(|i 9 | 2 ) were evaluated. 
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Parameter 


Value (LJ units) 


Value (SI units) 


h c 


6.2 ±0.4 e 


2.2 ■ 10~ 20 J 


(/to 


0.15 ±0.09 a~ 2 


0.42 am" 2 


k A /tl 


1.3 ±0.3 e/af 


3.6 • 10~ 20 J/nm 4 


CO 


-0.05 ± 0.02 a" 1 


-0.08 nm" 1 


kG 


-0.26 [-2.8-0] e 


-1-0-10" 20 J 


kx 


1.5 ± 1 e/af 


4.2 ■ 10~ 20 J/nm 4 


7a 


0.007 ± 0.01 e/a 2 


0.7 ■ 10 -22 J/nm 2 



TABLE I: Elastic constants of our model membrane as ob- 
tained from a fit of the fluctuation spectra of pure membranes 
to the elastic theory. Values in SI units are estimates based 
on the identification a t ~ 6A and e ~ 0.36 ■ 
for explanation). 



The results are shown in Fig. [2] The data illustrates 
the characteristic features of the spectra: The height fluc- 
tuations are Goldstone modes, hence the height spectrum 
diverges for small wavelength modes. In contrast, the 
monolayer thickness spectrum is limited by the equilib- 
rium thickness and tends towards a constant value for 
small wavevector values 70]. It exhibits a characteristic 
peak at q 2 a 2 ~ 0.4, corresponding to a soft peristaltic 
mode with wavelength ~ 10<7f. At small q 2 , the fluc- 
tuation spectra are dominated by bending deformations. 
For larger values of q 2 , the spectra are dominated by the 
protrusion modes and are equal for the height and the 
monolayer thickness fluctuations. 

The solid line in Fig. [2] (main frame) shows the fit to 
the elastic theory, Eqs. Q13p and p4[) . with fit parame- 
ters k c , CAo, kA/t^, k\, and 7a- The results of the fit are 
given in Table HI The elastic theory describes the data 
in an excellent way. This confirms the validity of the un- 
derlying assumptions, most notably, the decoupling ap- 
proximations (see above). The strength of the coupling 
between bending modes and protrusion modes can be es- 
timated by looking at the coupling parameter ^\/k\k c . 
Our fit gives j 2 /k\k c = 5 • 10~ 6 , which is indeed much 
smaller than unity. The inset of Fig. [2] shows the fit 
of the monolayer thickness spectrum to the Landau-de 
Gennes theory (solid line). To make the analyses com- 
parable, we have included a protrusion contribution and 
fitted the monolayer thickness deformations to 



(l^)l 2 ) 



k R T 



knT 



A{a + cq 2 ) 2{k x + lx q 2 )' 



(27) 



(cf. Eq. (fl~2"|) ) and the bending deformations to Eq. (fT3|) 
simultaneously, with fit parameters a, c, k c , k\, and 7a . 
Not suprisingly, the Landau-de Gennes theory cannot re- 
produce the peak at nonzero q in (\t(q)\ 2 ), hence it misses 
one important characteristic of the thickness spectrum. 

Turning back to the elastic theory, the analysis of fluc- 
tuation spectra yields the values of three elastic param- 
eters that are needed in the subsequent analysis: the 
bending rigidity k c , the compressibility modulus Hai and 
the extrapolated curvature £. The two remaining elas- 



tic parameters in Eq. (|15p are the spontaneous curvature 
Co, and the Gaussian rigidity kQ. Under the (admit- 
tedly bold) assumption that the two monolayer slabs can 
be treated as elastic continua subject to internal stress, 
these quantities can be calculated from the first and the 
second moment of the surface tension profile across the 
monolayers [73j via 



(28) 



k c c = - / dz ji n t(z) (z - Zq) 



k G = 2 dz 7int (z) (z-z ) 2 . (29) 
Jo 

Here 7i n t(z) denotes the intrinsic surface tension profile, 
which is defined as the difference between the normal and 
the tangential components of the local pressure tensor. 
The integration starts at the bilayer midplane z = 0, 
and the reference plane z — zq is the "inextensibility 
plane", i.e., the plane in which an infinitesimal volume 
element is not compressed or extended if the monolayer 
is bent. For symmetric bilayers with overall tensionless 
monolayers, one has J °° dz 7i n t(-z) = 0, and the result 
for Co is independent of zq. The value obtained for kc, 
however, depends sensitively on the choice of zq. 

In the simulation, the pressure tensor is obtained using 
the virial theorem, 



T> NkBT * _i_ 1 

Vap = ~^S aP + — 



(30) 



where rj is the position of particle i, Fj the force act- 
ing on this particle, N the number of particles, T the 
temperature, and V the volume. Because of the periodic 
boundaries, care must be taken to use a version of the ex- 
pression (f3"0")) that does not depend on absolute positions, 
e.g., by rewriting the contribution of pairwise forces as 
J2i<j( r f ~ r< j)Fij> etc - The pressure tensor of the whole 
equilibrated system is diagonal, T a p = SapP, where P 
is the applied pressure. Nevertheless, it varies spatially 
on the molecular scale. In order to measure the local 
pressure profile, we divide the system along the z-axis in 
slices of length dz = 0.125o\ The contributions of pair- 
wise forces to the total virial are parcelled out on the 
bins according to the Irving-Kirkwood convention [74| . 
i.e., they are distributed evenly on the line connecting 
the two interaction partners. The contributions of our 
only multibody interactions, the bond-angle potentials 
(Eq. (0|), are distributed evenly on the two participat- 
ing bonds. Alternatively, Goetz and Lipowsky [75| have 
proposed a method where the contribution of multibody 
forces is spread evenly on all lines connecting the par- 
ticipating partners (the bond-angle virials would then be 
distributed on triangles). We have also implemented this 
definition, and found that the effect on the pressure pro- 
file was negligible. The interfacial tension profile is given 
by 



l{z) = V zz {z)--{V xx {z)+V vy {z)). 



(31) 
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FIG. 3: Surface tension profile 7(2) for four different system 
sizes with N = 200, 288, 800, and 3200 lipids (top panel). 
The bottom panel shows the corresponding density profiles 
of solvent, head and tail beads in the smallest system (200 
lipids) for comparison. 



Figure [3] shows the surface tension profiles of a pure 
model membrane for four system sizes. Qualitatively, 
it exhibits the same features as the stress profile ob- 
tained with the similarly simple coarse-grained lipid 
model of Goetz and Lipowsky [75(. Most notably, the 
surface tension features a negative peak in the bilayer 
midplane (at z = 0), indicating that the monolayers 
are strongly bound to each other, in agreement with 
atomistic and coarse-grained simulations of DPPC bilay- 

ers HlJIzilrzl- 

Due to the height fluctuations, j(z) depends on the 
lateral system size: All profiles are broadened in large 
systems. In contrast, the expressions (|28t and (|29|) are 
based on a hypothetical intrinsic surface tension profile 
7int( 2 0- If such an intrinsic profile exists, the actual ten- 
sion profile 7(2) should be given by the convolution of 
Tint (z) with the distribution of interface heights W(z'), 



7(as) = Jdz' W{z') liat {z-z'). 



(32) 



In this case, one easily verifies that the integral To = 
Jdz 7(2) is still zero for tensionless membranes, and 
the second moment I^ = Jdz z 2 7(2) does not de- 
pend on the shape of the function W{z) for symmet- 
ric tensionless bilayers with Jdz z 7(2) = 0. This pre- 
diction can be used to test the convolution hypothesis, 
Eq. (|32p . The integral To is close to zero in all sys- 
tems as expected (r of /e = -0.018, -0.003, -0.04, and 
0.02 for the system with N = 200,288,800, and 3200 
lipids, respectively). The values of the second moment 



are T 2 /e = 2.8,2.7,2.0, and 2.8. Hence T 2 does not de- 
pend on the system size, which confirms the existence of 
an intrinsic tension profile. 



We note that the integral for the Gaussian rigidity 
still depends on the system size, since the lower integra- 
tion bound is finite. Both Eqs. (|2"8"j) and (|29p are only 
applicable in sufficiently small systems. Therefore, we 
proceed by evaluating them in the smallest system with 
N = 200 lipids. For the monolayer curvature, we ob- 
tain a small negative value, k c Co = —0.3 ± O.le/ot. This 
may seem surprising, given the fact that the heads in 
our model are larger than the tail beads, and that the 
tails have to tilt in the low temperature phase (the 
phase) in order to accommodate this mismatch. Indeed, 
the spontaneous curvature is found to be positive in the 
gel state (78[. At higher temperature, the tails disor- 
der and occupy more membrane area, and Co decreases 
and changes sign as a result. Negative curvatures have 
also been found in more realistic models of DPPC bi- 
layers 54]. The calculated result for the Gaussian rigid- 
ity depends on the position of the "inextinsibility plane" 
zq, which is not known unambiguously. Depending on 
our choice of zq, we obtain kc values that range between 
±2.8e. Among these, the positive values can be excluded: 
Positive Gaussian rigidity would imply that the bilay- 
ers tend to assume saddle-shaped configurations, which 
destabilizes flat bilayer structures on large scales and 
promotes interconnected structures, i.e., cubic phases or 
sponge-like structures. Our model bilayers remained flat 
for all system sizes. Therefore, only negative values of 
kc are physically reasonable. In previous work 32|, zq 
was taken to be the "neutral" plane where "f(z) crosses 
zero, i.e., zq = 2.604 in our system (cf. Fig. ©). Using 
this value, we obtain k& = — 0.26e. 



To summarize this subsection, the Landau-de Gennes 
theory does not describe the fluctuations of our model 
membranes very well. It misses a soft peristaltic mode 
in the monolayer thickness spectrum, which is clearly 
present in the simulation data. In contrast, the elastic 
theory fits the data excellently. We have extracted val- 
ues for the elastic parameters that describe the bending 
deformations in monolayers from the fluctuation spectra 
and the surface tension profiles of bilayers. They are 
given in Table fl] Using the identifications a t ~ 6A and 
e ~ 0.36 • 10 _20 J (see section 2), we can relate our model 
to real lipid bilayers. The values of our elastic parameters 
in SI units (see Table [TJ have the same order of magni- 
tude than those obtained based on all-atom simulations 
of DPPC 0, k c ~ 4 • 10- 20 J, k A /t 2 ~ 1.1 • 10- 20 J/nm 4 , 
C/to ~ 0.18nm" 2 [13, and c ~ -0.04 - -0.05 nm" 1 H, 
and those based on experimental estimates [79j], k c ~ 
5 - 20 ■ 10- 20 J, k A /4 ~ 6 ■ 10- 20 J/nm 4 , c 0.04 



nm , and kc/k c 



-0. 



9 




Distance from Protein [o ( ] 



FIG. 4: Radial membrane thickness profiles in the vicinity 
of inclusions with different hydrophobic thickness L and hy- 
drophobicity parameter e pt as indicated. Filled symbols show 
data for inclusions with fixed orientation along the 2-axis, 
open symbols correspond to unconstrained inclusions. The 
solid and dashed lines are fits to the elastic theory (Eqs. (|20p 
with (|17[) . and (|18[) ') with fit parameters ta and cq. 
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FIG. 5: Membrane thickness profiles in the vicinity of an in- 
clusion with hydrophobic thickness L — 4a t ,6at, and 8<7t, and 
hydrophobicity parameter e p t = 6, compared with fits to the 
Landau-de Gennes theory (dashed lines), to the elastic the- 
ory with fixed Co = —0.05 /at and kc = — 0.26e (solid lines), 
and to the elastic theory with the additional fit parameter cb 
replacing Co (dotted line). The grey lines indicate the range 
of the fit at fixed Co if Co and ka are varied within the error 
given in Table U 



B. Deformation of a Bilayer by a Single Protein 

Next we investigate the deformation of a bilayer by a 
single inclusion. To this end, we consider the radial pro- 
files of the membrane thickness 20, which we define as the 
mean z-distance between the head positions in the upper 
and the lower monolayer. Two parameters are varied, 
the hydrophobicity parameter e pt (cf. Eq. ([6])) and the 
hydrophobic thickness L of the inclusion. The results for 
e p t = 2 and 6 are shown in Fig. [3] for proteins with freely 
fluctuating orientations and for proteins with orientation 
constrained to the the z-axis. The deformation profiles 
induced by proteins with fixed and free orientiations are 
identical within the error. This is due to the fact that 
proteins with free orientations were hardly tilted - the 
tilt angles were always smaller than 9 ~ 0.08. In the 
following, we shall only show data for proteins with fixed 
orientation. 

Looking at Fig. 2J we first observe that the mem- 
brane thickness profiles are not strictly monotonic, but 
exhibit a characteristic over- or undershoot at the dis- 



tance r ~ 6<7t from the protein. Such a weakly oscilla- 
tory behavior has also been observed in previous coarse- 
grained (orT ] and atomistic (80| simulations of protein- 
induced membrane deformations. In our case, the wave- 
length of the oscillation is roughly ~ 10<7t, hence it can 
be related to the soft peristaltic mode in the fluctuation 
spectrum. 

The second observation is that the protein hydropho- 
bicity parameter e pt must exceed a certain value in or- 
der to produce classical hydrophobic mismatch. If e pt is 
too small, the protein effectively repels the lipids, and 
the membrane thickness is reduced at the surface of the 
protein regardless of the value of L. The hydrophobic 
section of the protein pins the membrane thickness for 
hydrophobicity parameters larger than e pt ~ 4. This can 
be rationalized by noting that e pt ~ 4 is the critical value 
where touching the protein surface is about as favorable 
for tail beads, from an energetic point of view, as being 
immersed in the bulk: The maximal contact energy of a 
tail bead in contact with a plane of tail beads is 4e. 

We proceed by fitting the profiles to the prediction 
of the analytical theories. Since the hydrophobicity of 
the inclusion must be larger than e p t > 4 in order to 
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L [at] tfl dG [<Jt] t% [at] Co [a t *] 

4 -0.94 ± 0.02 -0.66 ± 0.02 -0.11 [-0.12 - -0.10] 

6 0.3 ±0.02 0.18 ±0.02 0.05 [-0.03 - 0.06] 

8 1.44 ±0.02 0.93 ±0.01 0.22 [0.15 - 0.26] 



TABLE II: Parameters i^ dG and ifj (monolayer deformation 
at the surface of the inclusion) and c~o (renormalized curva- 
ture) obtained from fitting the deformation profiles for large 
hydrophobicity parameters e p t — 5 and 6 to the Landau-de 
Gennes theory (tfl dG ) and to the elastic theory (f^ and do). 
The remaining fit parameter in the Landau-de Gennes fit is 
the decay length £ = 2.0±0.1<Tt. The uncertainty of cb results 
from the uncertainty of ka- 



cleave the membrane, we focus on the data for e pt = 6. 
Fig. [5J compares them with fits to the prediction of the 
Landau-de Gennes theory (Eq. (fit))) and the elastic the- 
ory (Eq. PD|) with the boundary conditions Eq. (118[) ). 
In the Landau-de Gennes case, the three curves were fit- 
ted simultaneously with one common fit parameter £ and 
three separate fit parameters tn = i^ dG . Not surprisingly, 
the Landau-de Genens fit cannot reproduce the oscilla- 
tory component of the profiles; otherwise, the fit is quite 
reasonable (Fig. [5J dashed line). 

The thick solid lines in Fig.[5Jshow the fits to the elastic 
theory with tn =: tfj as sole fit parameter. None of them 
is satisfactory. Hence the "pure" version of the elastic 
theory, which explains the profiles in terms of bulk mem- 
brane properties only, is not sufficient. In contrast, the 
data can be fitted very nicely if we release the constraint 
on the value of cq, i.e., replace the spontaneous curva- 
ture by a renormalized curvature do (thin dotted lines in 
Fig. [5J solid lines in Fig. 2]) . The resulting fit parameter 
values are essentially the same for e pt = 6 (Fig. [5j and 
£pi — 5 (data not shown) and given in Table [Til 

From the fit results for tn, we can infer the effective 
hydrophobic thickness L c ff = 2(fo + £r) of the inclusions. 
We note that the exact relation between tn and the model 
parameter L is not a priori clear, since the lipid-protein 
potential is smooth and varies on the length scale at- 
In all cases, the values for L c g are reasonably close to 
L, i.e., well within ler t . The fit parameters for do, how- 
ever, deviate considerably from the spontaneous curva- 
ture Co = —0.05 ± 0.02cr t " 1 (see Table [JJ and depend on 
L. This clearly demonstrates that the local structure of 
the lipids that surround the inclusion indeed contributes 
to the boundary conditions, as has been discussed in the 
theory section (Eq. |HJ) ). 

A similar effect has been observed by Brannigan and 
Brown in a different coarse-grained model in Ref. [32| . 
and could be explained satisfactorily by the effect of non- 
constant lipid volume. The volume per lipid in that 
model varies substantially over a range of r. In our 
model, the lipid volume (i.e., the lipid density inside 
the membrane) is almost constant throughout the sys- 
tem. Fig. [5] (upper left) shows profiles of the lipid bead 




Distance [a ( ] Distance [a,] 

FIG. 6: Radial profiles of various quantities as a function of 
the distance from the center of an inclusion with hydropho- 
bic thickness L — 4a t (black circles), L — 6a t (stars), and 
L — 8a t (white squares) at e pt = 6. Upper left: bead den- 
sity in the bilayer. Upper right: Nematic order parameter 
for bonds. Middle left: Heads per area in the monolayers. 
Middle right: Hydrophilic shielding parameter (see text and 
Ref. [U|). Lower left and lower right: Monolayer overlap, 
evaluated according to two different prescriptions (see text 
for definitions). 



density in the membrane for different hydrophobic thick- 
nesses L at e pt — 6. Here the lipid bead density p; is 
defined as the number of lipid beads per area divided 
by the membrane thickness, and it is directly related to 
the local lipid volume vi via vi = n/pi with the chain 
length n = 7. Apart from an enhancement directly at 
the protein surface, which reflects the attractive protein- 
lipid interaction, and a very shallow depletion zone there- 
after, the lipid density is nearly constant. Moreover, the 
curves for different hydrophobic thickness L are almost 
identical. Hence lipid volume effects contribute at most 
an L-independent constant to the renormalized curvature 
HMD- 
Fig. [6] also displays profiles of other candidate quanti- 
ties that might affect the membrane properties and renor- 
malize the curvature term at the surface. The upper 
right panel shows the nematic order parameter for sin- 
gle bonds, S z = (3(u z /u) 2 — l))/2, where u refers to 
the bond vectors. The profiles show that the nematic 
order is enhanced in the vicinity of the inclusion. This 
is partly related to the increase of lipid density at the 
inclusion, since both quantities are correlated. The de- 
tails of the nematic order profile, however, cannot be ex- 
plained by the density profile alone. The ordering effect is 
highest for positively mismatched inclusions. Suprisingly, 
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negatively mismatched inclusions induce higher order at 
the boundary than hydrophobically matching inclusions. 
This may seem counter-intuitive at first, but becomes 
plausible in view of the fact that the monolayers also 
overlap in the vicinity of negatively mismatched inclu- 
sions. At larger distances from the inclusion, the bond 
order decays monotonically for hydrophobically match- 
ing and positively mismatched inclusions, and exhibits a 
non-monotonic dip for negatively mismatched inclusions. 

The middle panel in Fig. [5] shows two quantities that 
are related to the shielding of the hydrophobic membrane 
interior from the solvent. Since shielding is achieved by 
the heads, the areal head density gives a direct measure 
of the effectiveness of shielding. At constant lipid volume, 
however, the areal head density depends directly on the 
monolayer thickness t(r) and thus does not qualify as in- 
dependent field 6q(r) in the elastic theory, i.e., it cannot 
contribute directly to Eq. (|2"4"|) . Instead, we must con- 
sider a rcnormalizcd shielding parameter, such as, e.g., 
the "hydrophilic shielding parameter" introduced by de 
Meyer et al. in Ref. 5l|. It is defined as the areal head 
density divided by the areal tail density, normalized such 
that it becomes one far from the inclusion. In the origi- 
nal version of the clastic theory, this parameter would be 
constant. Our data show that the areal head density is 
enhanced directly at the surface of the inclusion for all L 
- an indirect consequence of the attractive protein-tail in- 
teraction. The subsequent behavior depends on the type 
of mismatch: For positively mismatched inclusions, the 
areal head density goes up, for negatively mismatched in- 
clusions, it goes down (Fig. [SI middle left panel). The hy- 
drophilic shielding parameter (middle right panel) shows 
the same behavior at intermediate distances from the in- 
clusion: Overshielding for positively mismatched inclu- 
sions, undershielding for negatively mismatched inclu- 
sions. Close to the inclusion, the curves turn around, in 
qualitative agreement with the observations of de Meyer 
et al. in Ref. [5l|. The resulting integral J drSq(r)/q 
(the K2 term in Eq. (|2"4"|) ) is roughly zero for L = 6<Jt and 
8at and positive for L = Aat- 

The bottom panel in Fig.[S]shows profiles of the mono- 
layer overlap, which measures the amount of interdigi- 
tation between monolayers. It has been calculated fol- 
lowing two different conventions: The left graph shows 
a chain-related overlap parameter originally introduced 
by Kranenburg et al. in Ref. [81| : It is defined as 
O chain = (2(l z — to)/l z ), where l z is the z-component of the 
end to end vector of chains and to the monolayer thick- 
ness. Far from the inclusion, this parameter is negative, 
indicating that the two monolayers are well separated. 
Close to the inclusion, it becomes positive for negatively 
mismatched inclusions - the inclusions pull the lipids in- 
wards and enforce a certain amount of intcrdigitation. 
For hydrophobically matching inclusions and for posi- 
tively mismatched inclusions, it remains negative and 
roughly constant. The right graph shows a monomer- 
related overlap parameter which is defined as the overlap 
integral of the density profiles of the two monolayers, 



O bcad = / dz(/9 t u a p if *{z) - p l °J"(z)), where p tail denotes the 
density of tail beads. The curves for O bcad are qualita- 
tively similar to those for O chain , except that they are 
shifted to positive values - even far from inclusions, the 
monolayer densities have some overlap in the center of 
the bilayer (cf. Fig. [3]). 

To summarize this section, we find that the thickness 
deformation profiles around single inclusions can be fit- 
ted reasonably well with the exponential law predicted 
by the Landau-de Gennes theory. The fit with the pure 
version of the elastic theory is not good, but the elastic 
fit becomes much better and far superior to the Landau- 
de Gennes fit, if we allow for the possibility of curvature 
renormalization in the vicinity of the inclusion. We have 
identified a number of quantities which could conceivably 
contribute to this renormalization, i.e., they vary sub- 
stantially close to the inclusion and they show a sizeable 
L-dependence. However, we cannot pinpoint an obvious 
single culprit. According to Table HH the renormalized 
curvature c exhibits an almost perfect linear dependence 
on the hydrophobic thickness L of the inclusion: The re- 
lation do - c » 0.078/er t + 0.0825(L - 2t )/cr t 2 describes 
the data in Table [IT] within 4%. None of the quantities 
shown in Fig. [6] produces such a linear behavior in an ob- 
vious way. We conclude that we have either still missed 
the truly relevant quantity, or that the observed linear 
relation is accidental and results from an interplay of 
various curvature-renormalizing factors. ^From a physi- 
cal point of view, it seems likely that most or all of the 
quantities discussed above will affect the local membrane 
properties, and most notably, the spontaneous curvature: 
Monolayer overlap will favor negative spontaneous cur- 
vature, since the area per tail increases. Variations in 
the hydrophilic shielding parameter will change the lo- 
cal pressure profile and hence affect the local curvature. 
Chain order will favor positive spontaneous curvature, 
since the structure in the chain region resembles that in 
the gel state, where the spontaneous curvature is large 
and positive. 



C. Effective Interactions between Two Proteins 

Finally, we turn to studying the membrane-mediated 
interactions between inclusions. We study the dilute pro- 
tein limit, i.e., we do not consider effects from multibody 
interactions between several proteins. In order to deter- 
mine the potential of mean force, we have carried out 
simulations of membranes containing two inclusions and 
determined their radial pair distribution function. The 
effective potential w(r) is 



w(r) = -k B T\ng(r), 



(33) 



where ks is the Boltzmann constant and T the tem- 
perature. Since the function g{r) varies over several or- 
ders of magnitude, we have used umbrella sampling and 
reweighting to determine it accurately at all distances r. 
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FIG. 7: Potential of mean force between two inclusions with 
hydrophobic thickness L — 4a t (negative mismatch, top 
panel), L = 6<7t (no mismatch, middle panel), and L = Sat 
(positive mismatch, bottom panel) for different values of the 
hydrophobicity parameter e pt as indicated. The inset in the 
bottom panel shows the interactions generated outside of 
the membrane (solvent-mediated depletion interaction and di- 
rect interaction) for hydrophobically matched inclusions (solid 
line), and the additional contribution of solvent-induced in- 
teractions at L = 4a t (dashed line) and L = Sa t (dotted line) . 



The resulting potential curves are shown in Fig. [7] for 
several values of the hydrophobicity parameter e pt and 
the hydrophobic thickness L. The most striking feature 
of these profiles is their distinctly oscillatory shape. The 
oscillations have a wavelength of roughly la t in all sys- 
tems, which indicates that they are caused by lipid pack- 
ing in the vicinity of the inclusions. Only for the lowest 
value of the hydrophobicity parameter, e p t — 1, the os- 
cillatory structure disappears. Here, we recover qualita- 
tively the behavior observed in the simulations of purely 
repulsive inclusions by Sintes and Baumgartner [52j | and 
predicted by the corresponding molecular mean-field the- 
ories [13, [H, : The potential of mean force exhibits a 
minimum at close inclusion distances followed by a shal- 



low maximum. 

Beyond r > 3.5 cr t , all layering minima obey the gen- 
eral rule that they become more shallow with increasing 
protein separation r and/or decreasing hydrophobicity 
parameter e pt . The layering effects are most pronounced 
at high e p t, indicating that lipids pack more tightly if 
they move closer to the protein surface. In contrast, the 
first potential minimum at r ~ 3 er t features the opposite 
behavior: It deepens as e p t decreases and disappears for 
high e p t. This minimum results from a number of effects 
that are not directly related to layering: (i) The direct 
protein-protein interaction and the solvent-induced de- 
pletion interaction between the hydrophilic protein sec- 
tion located outside of the membrane: This contribution 
is roughly constant and can be calculated analytically 
(inset of Fig. [7]). (ii) A depletion-type interaction in- 
duced by the lipids: By pushing the proteins towards 
each other, they maximize their translational and con- 
formational entropy. This effect is strongest at low e p t, 
where the protein and the lipids effectively repel each 
other, (iii) A bridging interaction induced by the lipids: 
At higher e pt , the lipids gain from being in contact with 
the proteins. Therefore, they tend to squeeze themselves 
between the proteins, pushing them apart, and the height 
of the first minimum goes up. 

The competition between the depletion interaction and 
the lipid bridging effect accounts for the phenomenon re- 
ported earlier (Fig. [1} , that the preferred arrangement of 
inclusions in the membrane depends on e pt : Weakly hy- 
drophobic inclusions tend to touch each other, whereas 
strongly hydrophobic inclusions favor a larger distance 
where they are separated by one single lipid layer. 

The influence of hydrophobic mismatch on the poten- 
tial of mean force can be assessed by comparing the po- 
tential curves for different hydrophobic thickness L. On 
the one hand, hydrophobic mismatch affects the local 
features of the potential - the packing interaction (the 
strength of the layering) and the effective contact en- 
ergy of proteins. On the other hand, it also contributes 
a smooth attractive interaction for mismatched proteins 
with L = Act t and 8at, which superimposes the oscillatory 
packing interaction at r > 4<j t - It is worth noting that 
the sign of the additional term is independent of the type 
of mismatch - it is attractive for both positively and neg- 
atively mismatched inclusions. This smooth long-range 
contribution to the potential can now be compared with 
the predictions of the analytical theories. 

Fig. [5] shows the corresponding curves for the Landau- 
de Gennes theory (thin lines) and the elastic theory (thick 
lines). They were calculated numerically by minimizing 
the free energy - Eq. ([9]) or Eq. (fT5|) with cq replacing Co 
- for systems containing two inclusions at given distance 
r, with the boundary condition (f> = tR at the surface 
of the inclusion. The model parameters were taken from 
Table [J and [TTJ In the Landau-de Gennes calculation, 
the parameter 4a in Eq. §§§ was identified with the re- 
duced area compressibility coefficient fc^/io m Table U 
This may seem inconsistent, since the latter was origi- 
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FIG. 8: Effective interaction potential between two inclusions 
according to the Landau-de Gennes theory (thin lines) and 
the elastic theory (thick lines) for inclusions with different 
hydrophobic thickness L as indicated. The thick grey line 
shows the simulation data for L — 4<r t , e pt — 6 for comparison. 



nally determined from a fit of the fluctuation spectra to 
the elastic theory. Unfortunately, the fit of the spectra to 
the Landau-de Gennes theory (Fig. did not produce 
dependable parameters a and c. The value of ItA in Ta- 
ble |T] is compatible with independent simulation data on 
the lipid area increase at finite surface tension [82| and 
was thus considered to be more reliable. The value for c 
then follows from c = £ 2 a, using the value £ = 2.0<7t de- 
termined from the fit to the membrane thickness profiles 
around single inclusions. 

To calculate the free energy, we have discrctizccl the 
corresponding integrals in real space using a square grid 
with spatial discretization parameter h and a second or- 
der difference scheme to evaluate the derivatives. The 
boundary condition was implemented by setting (j> = tj{ 
inside the inclusion. (The other boundary condition 
of the elastic theory, Eq. (fT5|) . follows automatically 
from the energy minimization). The energy was mini- 
mized via a steepest descent method, using a relaxation 
scheme introduced in Ref (83[. The final accuracy was 
/ d 2 r\5F/6(l)\ < 1Q- 6 . The curves shown in Fig. [5] were 
obtained using the spatial discretization h = 0.25at and 
a system of size 30 x 20of with periodic boundary con- 
ditions, which corresponds to the situation in the Monte 
Carlo simulations. Calculations for h — 0.5<7t and system 
sizes up to 70 x 50of produced the same curves, hence 
discretization errors and finite size effects are negligible. 

Both the Landau-de Gennes theory and the elastic the- 
ory predict an attractive interaction at distances r < 6174, 
in agreement with the trend observed in the simulation 
data. For larger distances, the elastic theory predicts 
the potential of mean force to become positive and ex- 
hibit a peak at r ~ 8a t . The simulation data show no 
indication for the existence of such a positive peak. Ap- 
parently, the elastic theory fails to reproduce the trend of 
the simulation data at large distances, which is surpris- 



ing, since this is where one would expect it to work best. 
The potential predicted by the Landau-de Gennes the- 
ory is negative and attractive everywhere and thus better 
compatible with the data. However, this should not be 
overrated, given the overall poorer performance of this 
approach when looking at pure membranes and mem- 
brane interactions with single inclusions. The weakly os- 
cillatory behavior of the potential of mean force in the 
elastic theory is generated by the soft peristaltic mode in 
the fluctuation spectrum, which has been shown to leave 
a clear signature in the shape of the distortion profiles 
around single proteins. The simulation data suggest that 
the effect of this mode on the lipid-mediated interactions 
between two proteins is destroyed by some yet unknown 
mechanism. 



IV. SUMMARY AND CONCLUSION 

To summarize, we have determined protein-membrane 
interactions and calculated lipid-mediated interactions 
between proteins in a simple generic molecular model for 
lipid bilayers, and compared the results to the predictions 
of two popular continuum theories. Whereas the effect 
of the protein-membrane interactions on the deformation 
profiles (Fig. 0]) can be described very nicely by a theory 
that essentially treats the membrane as a pair of cou- 
pled elastic sheets (monolayers), the local lipid structure 
clearly dominates the shape of the membrane-induced 
protein-protein interactions. 

We note that the specific shape of these packing in- 
teractions depends sensitively on the microscopic details 
of the system. All liquids have a local liquid structure, 
and packing effects will clearly also be present in real 
membranes. However, their contribution to the poten- 
tial of mean force will differ from that in our model. In 
particular, the amplitude of oscillations will presumably 
be much smaller, since packing effects are most likely 
overestimated in our simple Lennard- Jones bead model. 
Hence the potentials of mean force shown in Fig. [7] can- 
not be related to the potentials of mean force of proteins 
in membranes in real systems in a quantitative sense. 
At a qualitative level, however, some conclusions can be 
drawn: 

We have identified several factors that may contribute 
to the effective potential between cylindrical proteins in 
our simple model membranes: (i) Direct protein-protein 
interactions, (ii) depletion interactions, (iii) lipid bridg- 
ing interactions, (iv) packing interactions, and (v) elas- 
tic interactions. The interplay and competition of these 
factors determine the final, most favorable arrangement 
of proteins in the bilayer. Their relative importance is 
determined by the hydrophobicity of the proteins and 
the hydrophobic mismatch. The most dramatic effect of 
hydrophobic mismatch in our system can be associated 
with its influence on the factors (ii)-(iv). Thus we do 
observe a "hydrophobic mismatch interaction" , but the 
dominant contribution to this interaction seems to re- 
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lated to local reordering effects, and not to the elasticity 
of the monolayers. Intriguingly, the elastic theory does 
not even describe correctly the smooth long-range part 
of the interaction. 

As a general trend, the interaction tends to be attrac- 
tive i.e., it is most favorable for proteins to cluster. This 
is trivially the case for purely oscillatory interactions, but 
it also seems to hold for the additional smooth contribu- 
tion, regardless of the type of mismatch. The observa- 
tion is consistent with the findings of de Meyer et al. [5l| 
- the potentials of mean force presented in that paper 
are also always attractive, except for proteins with very 
large diameters. Whereas protein clustering may some- 
times be desirable, it also has negative effects - for exam- 
ple, it reduces the mobility of the proteins in the mem- 
brane, it makes them less accessible for other proteins, 
etc. Our results thus raise the question whether the ef- 
fects reported here have to be counteracted by external, 



not membrane-related protein interactions, or whether it 
is possible to identify mechanisms that induce repulsive 
membrane-mediated interactions between proteins, and 
stabilize protein dispersions, e.g., in mixed multicompo- 
ncnt membranes. 
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